The Schottky defect formation energy in MgO calculated by diffusion Monte Carlo 
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The energetics of point defects in oxide materials plays a major role in determining their high- 
temperature properties, but experimental measurements are difficult, and calculations based on 
density functional theory (DFT) are not necessarily reliable. We report quantum Monte Carlo 
(QMC) calculations of the formation energy Es of Schottky defects in MgO, which demonstrate 
the feasibility of using this approach to overcome the deficiencies of DFT. In order to investigate 
system-size errors, we also report DFT calculations of Es on repeating cells of up to ~ 1 000 atoms, 
which indicate that QMC calculations on systems of only 54 atoms should yield high precision. The 
DFT calculations also provide the relaxed structures used in the variational and diffusion Monte 
Carlo calculations. For MgO, we find Es to be in close agreement with results from DFT and from 
model interaction potentials, and consistent with the scattered experimental values. The prospects 
for applying the same approach to transition metal oxides such as FeO are indicated. 
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The quantum Monte Carlo (QMC) technique [1] is im- 
portant in condensed matter science, because it is gen- 
erally much more accurate than density functional the- 
ory (DFT) and allows one to make accurate predictions 
for problems where DFT fails [2-6]. ft is usually com- 
petitive in accuracy with high-level quantum chemistry 
methods, but it has the advantage of being practicable for 
large systems containing hundreds of atoms. Recently, ef- 
forts have been made to apply QMC to oxide systems, 
including transition metal oxides [7,8], for which DFT 
gives poor predictions for magnetic properties, phonon 
frequencies and other properties [9] . We have recently re- 
ported [10] a QMC study of bulk MgO for which we find 
excellent agreement with experiments for the bulk lattice 
parameter and bulk modulus, provided appropriate cor- 
rections are made. We report here what we believe to be 
the first QMC calculation of an oxide lattice defect en- 
ergy, namely the Schottky formation energy Es in MgO. 
Defect energies in oxide materials are technologically im- 
portant for applications ranging from high-temperature 
superconductors to radioactive waste disposal, but are 
also very difficult to measure experimentally. We shall 
show that our results for Eg in MgO are consistent with 
available experimental data, as well as supporting earlier 
DFT predictions [11]. 

The Schottky energy Es is the energy required to form 
a cation and anion vacancy pair, and governs the thermal 
equilibrium concentration of vacancies [12]. The calcu- 
lation of Es in ionic materials has a very long history, 
going back to the very early work of Mott and Little- 
ton [13-15]. DFT calculations on defect formation and 
migration energies in oxides first became possible in the 
early 1990's, and are now routinely performed. However, 
particularly in transition metal oxides, the reliability of 
DFT calculations is questionable. 



QMC calculations are usually performed in two 
stages [1]. In the first, known as variational Monte Carlo 
(VMC), a trial many-body wavefunction is constructed 
as a product of a Slater determinant of single-electron 
orbitals and the so-called Jastrow factor which explicitly 
accounts for electronic correlation. Since VMC by itself 
is not usually accurate enough, the second stage is to use 
the many electron wavefunction produced by VMC in dif- 
fusion Monte Carlo (DMC), which improves the ground 
state estimate by performing an evolution in imaginary 
time. This would yield the exact ground-state energy but 
for the fact that the trial wavefunction fixes the nodes of 
the many-body wavefunction. Because of this, the DMC 
energy is an upper bound to the true ground-state energy, 
but for systems with a large band gap the difference is 
expected to be very small. 

QMC calculations of defect energies in any material 
are challenging, and their feasibility is not obvious, for 
several reasons. First, it is not yet routinely possible to 
perform structural relaxation with QMC, and relaxation 
around defects produces a very large energy lowering in 
ionic materials [13-15]. Second, defect energies must be 
obtained as a difference of two large energies, both of 
which suffer the statistical errors inevitable with Monte 
Carlo methods. Third, there is a system size error associ- 
ated with the limited size of the periodically repeated cell 
used in condensed-matter QMC calculations, and the dif- 
ficulty of going to very large cell sizes makes it difficult to 
assess the error. To overcome these problems, we employ 
DFT calculations in tandem with QMC. In particular, 
the relaxed ionic positions used for QMC calculations on 
the defective crystals are taken from DFT calculations. 

The overall strategy for our calculations is as follows. 
The Schottky energy Es is defined to be the energy 
change when a Mg 2+ ion and an O 2- ion are removed 
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from the MgO perfect crystal, the resulting defective 
crystal is allowed to relax, and the removed ions are re- 
placed to form new perfect crystal. The two vacancies 
formed in this process are supposed to be very far apart, 
so that there is no interaction between them. In calcu- 
lating Eg, either by DFT or by QMC, it is preferable 
to perform calculations in which no more than a single 
vacancy, either a Mg 2+ vacancy or an 2 ~ vacancy, is 
present. We avoid doing calculations with both vacan- 
cies present in the same system, because with the cell 
sizes that can be achieved in practice, such calculations 
would suffer from a large unwanted interaction between 
the vacancies. Denoting by Em{v + , v~) the relaxed total 
energy of a crystal containing N cation lattice sites and N 
anion lattice sites, and with v + cation vacancies and v~ 
anion vacancies, we therefore express the Schottky energy 
as E s = E N (l, 0) + E N (0, 1) - [2(JV - l)/N] E N (0, 0). ' 

The energies En(1,0) and £^(0,1) both refer to pe- 
riodic systems in which the supercell has a net charge. 
To make these energies well defined, we must assume 
that the systems have been rendered electrically neu- 
tral by the introduction of a uniform background charge, 
as described in earlier papers [11,16,17]. The effect of 
the background charge densities is to make the zero- 
wavevector terms in the Coulomb energy finite rather 
than infinite, but these charge densities do not appear 
explicitly in any other way in cither DFT or QMC. 

In this way of doing the calculations, the interaction 
of the charged vacancies with their periodic images gives 
finite-size corrections AE to the total energy, which scale 
as L , where the length L characterises the dimensions 
of the supercell. As is well known [11,16], these correc- 
tions can be quite accurately approximated by the for- 
mula AE ~ aq 2 /eoL, where eo is the static dielectric 
constant of the material, q is the net charge of the va- 
cancy, and a is the appropriate Madelung constant. Us- 
ing this formula, we can subtract off the leading finite-size 
corrections, and greatly accelerate the convergence with 
respect to supercell size. In applying this correction pro- 
cedure in the present work, we take e for bulk MgO from 
experiment [20]. 

The present DFT calculations are performed using the 
VASP code [21] on a wide range of periodic cell sizes, 
as described above. These calculations allow us to de- 
termine the cell size needed to obtain converged results. 
The calculations employ the projector augmented wave 
method [22,23] for the interactions between the valence 
electrons and the ions, and the local density approxima- 
tion for electronic exchange and correlation. In all the 
DFT calculations, the entire system is fully relaxed, so 
that the maximum force on any ionic core is less than 
4 x 1(T 4 cV A" 1 . 

Detailed descriptions of VMC and DMC and of the 
CASINO code used for all the QMC calculations have been 
reported elsewhere [1,24]. Our trial wavefunctions 
have the usual Slater- Jastrow form \&t = D^D^ exp(J), 



where and D- 1 are Slater determinants of up- and 
down-spin single-electron orbitals, and exp( J) is the Jas- 
trow factor describing correlations between electrons. 
The function J is a sum of parameterised one- and two- 
body terms, the latter being designed to satisfy the cusp 
conditions. The free parameters in J are determined by 
requiring that the variance of the local energy in VMC be 
as small as possible. The many-body wavefunction rep- 
resents explicitly only valence electrons, whose interac- 
tions with the ionic cores are described by pseudopoten- 
tials. The Hartree-Fock pseudopotentials used here are 
the same as those used in our recent QMC work on bulk 
MgO [10]. The single particle orbitals have been taken 
from DFT-LDA calculations with the same pseudopo- 
tentials using the pwscf code [18]. The basis set used 
for the representation of the single-particle orbitals in 
QMC is the recently described B-spline or blip-function 
basis [25], which is closely related to plane waves, but 
is computationally far more efficient for QMC. Since we 
cannot perform structural relaxations with QMC, the re- 
laxed ionic positions used in our QMC calculations on the 
defective systems are taken from our DFT calculations. 

In order to suppress statistical bias, QMC calculations 
need to be run with a large population of "walkers" , and 
this makes it efficient to run them on large parallel ma- 
chines. The present calculations were performed on the 
HPCx machine at Daresbury using between 128 and 320 
processors, with a target number of 640 walkers. 

In Fig. 1 we display the value of Schottky energy cal- 
culated using DFT on various cell sizes, going up to 1024 
atoms [19]. The results include the Coulomb correction 
mentioned above [16]. Two sets of calculations are re- 
ported in the figure: one performed by sampling the 
Brillouin Zone (BZ) at the T-point only, and the second 
using a 2 x 2 x 2 Monkhorst-Pack [26] grid. The two sets 
of results become essentially indistinguishable for cells 
containing 128 atoms or more and converge very quickly 
to the value of 6.76 eV. The error in the Schottky en- 
ergy obtained from calculations on cells containing only 
54 atoms is somewhat less than 0.2 eV, but we note that 
cancellation of BZ errors makes the results obtained with 
54 atoms and T-point only sampling already converged 
to within ~ 0.07 eV. 

DMC calculations have been performed on cells con- 
taining 54 atoms, using a time step of 0.005 a.u. for over 
50,000 steps. With this length of simulations, total en- 
ergies for the perfect MgO crystal and the crystals with 
one Mg 2+ or one 2 ~ vacancy were obtained with sta- 
tistical errors of - 0.23 eV, - 0.32 eV and - 0.27 cV 
respectively. The value of the Schottky energy with the 
error bar obtained by combining the errors on the var- 
ious components is 7.5 ± 0.53 eV. For completeness, we 
also report the value of the Schottky energy of 6.99 eV 
obtained using DFT-LDA with the same pseudopoten- 
tials used in the DMC calculations. Although there are 
accurate experimental data for the migration energies of 
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cation and anion vacancies in MgO, the Schottky energy 
itself is experimentally uncertain, with measured values 
spanning the range 4 — 7 eV [27]. Our QMC value is con- 
sistent with this, and also with earlier predictions for Eg 
from both DFT [11] and calculations based on interac- 
tion models [28,29], all of which give values in the range 
6.5 - 7.5 eV. 

These results demonstrate the technical feasibility of 
using high precision QMC calculations to study the en- 
ergetics of defects in oxide materials. This is encouraging, 
because it suggests the possibility of using QMC to calcu- 
late the formation, association and migration energies of 
other kinds of defects, including impurities in oxide mate- 
rials. In the present work, the QMC result for the Schot- 
tky energy agrees with the DFT predictions to within the 
QMC statistical error, but this is expected because DFT 
is known to give a good description of most properties of 
MgO. However, this will not be true of strongly correlated 
materials such as transition metal oxides, for which DFT 
predictions are often poor. We are currently attempting 
to extend our QMC calculations to the important oxide 
FeO. 

A technical difficulty that is clear from the present 
work is that very long QMC runs are needed to reduce 
the statistical error on the defect energy to an acceptable 
level. This difficulty can be mitigated by improving the 
quality of the trial wavefunction. In the present case, we 
might do this by allowing the Jastrow factor to be differ- 
ent in the vacancy region. The use of recently developed 
methods for improving the scaling of the computer effort 
with system size [30-32] (so-called O(N) QMC) will also 
help in future work. But the fundamental problem is that 
there is no cancellation of statistical noise on the large 
energies that are subtracted. Correlated sampling may 
perhaps help with this, but the solution may lie also in 
the close connection between O(N) and embedding [33] 
that has already explored within tight-binding and DFT 
methods. 
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FIG. 1. Schottky vacancy formation energy (eV) in 
MgO calculated in the local density approximation of 
DFT for repeating cells containing different number of 
atoms. Filled circles connected by solid lines (filled squares 
conected by dotted lines) show results obtained with T-point 
(2x2x2 Monkhorst-Pack) Brillouin zone sampling. 



4 



